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tJQ; ABSTRACT 

We present a method for performing data-driven simulations of solar active region 
formation and evolution. The approach is based on magnetofriction, which evolves 
the induction equation assuming the plasma velocity is proportional to the Lorentz 
force. The simulations of active region coronal field are driven by temporal sequences 
of photospheric magnetograms from the Helioseismic Magnetic Imager (HMI) instru- 
■ ment onboard the Solar Dynamics Observatory (SDO). Under certain conditions, the 

data-driven simulations produce flux ropes that are ejected from the modeled active 
region due to loss of equilibrium. Following the ejection of flux ropes, we find an en- 
hancement of the photospheric horizontal field near the polarity inversion line. We also 
present a method for the synthesis of mock coronal images based on a proxy emissivity 
calculated from the current density distribution in the model. This method yields mock 
coronal images that are somewhat reminiscent of images of active regions taken by in- 
struments such as SDO's Atmospheric Imaging Assembly (AIA) at extreme ultraviolet 
wavelengths. 



1. Introduction 

Understanding the structure and evolution of the solar coronal magnetic field is an important 
component in understanding space weather. These dynamics are a consequence of the transport 
of energy through the photosphere and into the chromosphere and corona, and of the subsequent 
reorganization by the coronal magnetic field in response to this energy input. Some of this energy 
is used to accelerate the solar wind, some is used to heat plasma trapped in closed coronal loops 
(after which it is emitted as radiation), and some is used to power flares and eruptive phenomena 
as coronal mass ejections (CMEs). 



Coronae above act ive regions have received particular attention (e.g., iRegnier &: Priestl 12007; 



Kazachenko et al.ll2012l ) because these are the sites that produce the most numerous and strongest 
flares, and are most often associated with CMEs and particle acceleration events. Additionally, 
active region coronae are copious emitters in the extreme ultraviolet (EUV) and x-ray wavelengths 
of light, with emissions typically taking the form of loop-shaped structures that are assumed to 
trace out the trajectories of coronal magnetic field lines. 
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Despite considerable progress in observing the structure and evolution of the solar corona, the 
root causes of many phenomena remain elusive. Much of this uncertainty stems from the fact that 
it is difficult to accurately map out the three-dimensional geometry of the coronal magnetic field, 
and to observe its temporal evo lution. Direct measurements of the coronal field are somet i mes 
avail able off the solar limb (e.g., iLin et al.1 12004 ; iBrosius fc White! 120061 ; iTomczvk et alj|2008l : iLm 
2009), but these measurements often lack sufficient spatial or temporal resolution, and suffer from 
line-of-sight confusion, thereby making proper interpretation somewhat challenging. 

As a result of these difficulties, much effort has been put toward modeling the coronal magnetic 
field using photospheric magnetograms and/or coronal loop imagery as constraints. Many models 
make use of photospheric measurements of the magnetic field as boundary conditions. These 
models include potential (current-free) field extrapolations which provide a general idea of the 
global connectivity, force-free models which accommodate currents but assume the corona is in 
static force balance, and more physically realistic magnetohydro dynamic (MHD) models which 
capture many more of the important dynamical processes. 

Of these various coronal magnetic field models, only the potential-field source-surface (PFSS) 
model can presently be computed in real time in a forecasting (predictive) capacity. However, 
because the active region coronae of particular interest are those that likely contain significant 
currents, the applicability of the PFSS model for understanding energy buildup and release in 
active regions is limited, restricting the PFSS model to instances where only visualizing the large- 
scale (global) geometry of the coronal field, or p roviding contextu al magnetic fields surrounding 
active regions during relatively quiescent periods (jRiley et al.ll2006l ). is needed. 



At the other end of the spectrum of coronal field models are MHD simulations (e.g. jMikic et al 



1999). These simulations are the most sophisticated in terms of physical realism, but, even with 
today's supercomputing technology, remain computationally expensive for active-region-sized do- 
mains at appreciable resolution. Additionally, because not all of the necessary photospheric bound- 
ary conditions (e.g., gas pressures and electric fields) are measured directly, and because they 
still require parameterized models to bypass treatment of the microphysics (e.g., coronal heating 
functions) which are not always well known, even MHD models are not free from assumptions. 

In between the PFSS and MHD models lie the class of nonlinear force-free field (NLFFF) 
models of the corona. NLFFF models contain no dynamics, but do allow for field- aligned electric 
currents and thus enable measurements of magnetic free energies and magnetic helicity. Fast al- 
gorithms that extrapolate the magnetic field upward into the corona have been developed, though 
with mixed results in terms of their ability to r eliably reproduce the observed coronal field struc- 



tures (ISchrijver et al 



2008 



DeRosa et al. 



2009). Yet, because the electric current systems that 



occupy the corona tend to be isolated in and around active regions, NLFFF models seem particu- 
larly well suited for detailed studies of the build-up and release of energy within the corona. 



One NLFFF solution technique, the magnetofrictional (hereafter MF) scheme (jYang et al 



1986 



Craig &: Sneydlll986l ). involves integrating a model forward in such a way as to reduce the 
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magnetic stress in the model. This is achieved by assuming that the inductive velocity is parallel 
to the Lorentz force, which over time causes the magnetic field in the model to relax to a force-free 
state. This method has previously been employed as a way to create a NLFFF from a non- force- 



domain ( 


Valor i et al. 


2005, 


2007) or via the insertion of a flux-rope structure into an otherwise 


potential field (e.g., 1 


ran Balleeooiien 2004: Bobra et al. 2008: Savcheva & van Balleeooiien 2009: 


Savcheva et al. 


2012). 



Alternatively, one may drive a MF model by allowing the lower-boundary data to change with 

time. Long-term studies of the formation and subsequent evolution of fila ment channels using MF 

schem es have been performed bv lvan Ballegooiien. Priest. Mackavl (|2000l ) and lMackav. Gaizauskas. &: van Bailee 
(|2000l ). in which it was found that potential-field models did not result in sheared arcades asso- 
ciated with filament channels. Instead, evolving a MF model over multiple rotations was needed 
to recover the skew of coronal fields observed above the polarity inversion lines along which fila- 
ments appear. Agreement between observations of filament handedness and their modeled coun- 
terpa rts was found to improve as the model was run for longer (multiple months) periods of 



time (jYeates. Mackav. van Ballegooijen 



2007 



2008 



Yeates fc Mackavil2009bl : lYeates. Constable, Martens 



2010 ) . indicating that the build-up and transport of energy over these time scales is likely an im- 
portant factor. The same evolving MF s cheme was used in studies investi gat ing the processes 
and frequencies by w hich filaments lift off (jMackav &: van Ballegooii er] l2006al ]r] |; lYeates Mackavl 



2009a 



Su et al. 



2011 



the results of which showing that the changing magnetic geometry associ- 
ated with neighboring bipolar active regions is critical in determining how the upward ejection of 
flux ropes occurs. 

Until now, detailed NLFFF models of active region coronae have primarily been constructed 
only from single, isolated vector magnetograms, and therefore incorporate no information about 
the prior evolution of the photospheri c and coronal magnet ic field. With the advent of the He- 
lioseismic and Magnetic Imager (HMI; IScherrer et ai1l2012al lbl) on board NASA's Solar Dynamics 
Observatory (SDO), time series of photospheric vector magnetogram data are now available at 
relatively high resolution (0.5" pixels) and cadence (every 12 min), with unprecedented spatial and 
temporal coverage. Such data enable detailed models of active regions to be constructed, wherein 
the MF scheme is utilized to advance the models in time in response to driving from time se- 
ries of photospheric magnetograms. This technique provides a way to accommodate the (possibly 
widespread) instances where the prior history of the state of the coronal magnetic field factors into 
the determination of its future state. 

We show in a series of articles that it is feasible to use a MF approach with a time-evolving 
lower boundary condition to model the coronal evolution of an active region over a week-long period 
of time, and that ejections of magnetic flux can be driven in a model corona using this approach. 
In this article, the first in the series, we outline the numerical methods and demonstrate the scheme 
on NOAA Active Region (AR) 11158 using idealized boundary conditions. AR 11158, which was 
on disk in mid Feb. 2011 and produced several major flares, is followed for a significant fraction of 
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its disk passage. 



2. Method 



Our ap proach is based on th e magnetofriction method, which was introduced by lYang et al 



(jl986l ) and ICraig fc Sneydl (|1986l ) to examine the relaxation of magnetic configurations toward 
force-free states. All MF models share a common assumption, namely 



-j x B. 

v 



(1) 



In other words, the plasma velocity v is everywhere proportional to the Lorentz force F = (47r) _1 j x 
B, where the current density is defined here as j = V x B, and v is a specified magnetofrictional co- 
efficient. Given this chosen form for the velocity field, the magnetic field B is evolved in accordance 
with the magnetohydrodynamic induction equation 



8B_ 

dt 



V x (vxB — rjj), 



(2) 



where r\ is the magnetic resistivity. 



At first glance, the assumption represented by Eq. ([T]) may appear ad hoc and somewhat 
unphysical. However, there are a number of properties which make this assumption desirable for 
modeling the evolution of magnetic fields in a magnetically dominated (low plasma f3) regime. First 
of all, the choice of a velocity parallel to the Lorentz force means the magnetic energy in the volume 
of interest evolves according to 



-(j x By + r)j 



72 



(3) 



where S = — (y x B — rjj) x B is the Poynting flux of magnetic energy. Since the right-hand side 
of this equation is never positive (both v and rj > 0), it is easy to see that, in the absence of a net 
Poynting flux (i.e. no energy injection), the total magnetic energy (M = J y B 2 /8irdV) inside the 
volume must decrease monotonically or stay constant in time. 

By itself, the property of M < is desirable but not necessarily meaningful for studying the 
temporal evolution of coronal magnetic fields. Another property of the MF method is that, by 
virtue of the fact that the induction equation is used to advance B in time, the topology of the 
field is preserved (except in the case of high magnetic diffusion). In order w ords, under ideal MF 
evolution, important to pological quantities such as relative magnetic helicity (jBerger fc Fieldlll984l ; 
Finn &: Antonsenlll985l ) remain conserved. In contrast, this property does not necessarily hold for 
successive magne tic field configura tions generated from NLFF F extrapolation schemes such as the 



Grad-Rubin (e.g. I Wheatland! 120061 ) and the optimization (e.g. IWheatland et al.ll2000l : IWiegelmann 



2004J) methods. Force-free field extrapolations based on these two approaches treat the construction 



- 5 - 



of model fields at two times t\ and ti as c ompletely in d ependent p rocesses. The application of MF 
to the construction of force- field fields by IValori et al.l (|2005l . 120071 ) is similar in that the force- free 
field obtained at time t\ does not influence the solution at time ti- In this regard, our application of 
MF to modeling the time evolution of coronal fields should be clearly distinguished from previous 
efforts that apply MF for force-free field extrapolation. 

One way to physically interpret the MF scheme is to think of Eq. ([1]) as a simplified momentum 
equation which specifies the velocity without regard to inertial effects. Another way to interpret 
the MF scheme is as having a non-linear diffusion term for the magnetic field. We note that 
magnetic field evolution by MF has the same mathematical character as ambipolar diffusion. Both 
effects contr ibute to an electric field in the induction equation proportional to B x (j x B). As 
reported by iBrandenburg fc Zweibel ( 19 94). ambi polar diffusion leads to the formation of sharp 
current sheets (see also ICheung &; Cameron! I2012I ). What is important to note, however, is that 
both ambipolar diffusion and magnetofriction do not permit reconnection (their electric fields are 
perpendicular to B). In order for the topology to change, Ohmic diffusion (which is present in our 
model) needs to operate in locations of current layers. 



2.1. Implementation 



Our implementation of the MF scheme follows closely the implementation bv lvan Ballegooiien et al 



(|2000l ). who used the method to study the formation of filament channels. The major difference 
between the two methods is in the implementation of the bottom boundary condition, which will 
be discussed in section 12.3.21 

As an alternative to solving for the magnetic field B in time, the method solves the induction 
equation for the vector potential A, 

dA 

— =iJ x B-r]j, (4) 

where by definition B = V x A. We use a staggered mesh such that Cartesian components of the 
vector fields A, B and j (here defined as j = V x B) are defined in the following way: A x and j x 
are defined at the midpoints of ce ll edges p arallel to x and B x is defined at the centers of cell faces 



with normal vectors parallel to x (|Yedll966l ). The same arrangement applies to the other Cartesian 



directions. Numerical values of all components of v are defined at the cell vertices. 

Let superscripts i G [0, iV^], j € [0, iVy] and k € [0, N z ] denote the coordinate locations of cell 
vertices in the x, y and z directions, respectively. The length of the cell edges in the three directions 
are Ax, Ay and Az. In order to evaluate B = V x A, second-order spatial derivatives of A of the 
form 
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Az 

are used. Evaluation of j also employs second-order spatial derivatives: 



(6) 



i+Uk B z -B z 
Jx ~ Ay (7j 

Az ' (8) 

Derivatives for the y and z directions use the same scheme. To evaluate the magnetofrictional ve- 
locity v at cell vertices, linear interpolation of j and B is used. Furthermore, the same interpolation 
scheme is used to evaluate the v x B term in Eq.(j4|). 

For the simulations presented here, the functional form of the frictional coefficient is given by 

v = v B 2 (l-e-'' L )- 1 , (9) 

where vq = 10 s Mm -2 , L = 10 Mm a nd z is the height above the bottom boundary (i.e. photo- 
sphere) . As pointed out by iLowl (|2010l ) , some past models t hat use the magnetof riction approach 
have adopted a constant magnetofrictional coefficient v (e.g. ICraig h Snevdll2005l ) while assuming 
a line-tied, rigid bottom boundary where flows vanish. This could lead to an undesirable mismatch 
between the boundary condition and the magnetofrictional velocities that in turn result in unphys- 
ical current sheets near the bottom boundary. The height dependent profile of v in Eq.([9]) is chosen 
such that magnetofrictional velocities smoothly transition to zero towards z = 0. 

It is worthwhile to estimate the relevant timescales based on the above choice of the mag- 
netofrictional coefficient. Since the MF induction equation is nonlinear, it is difficult to estimate 
relaxation timescales with general validity. To estimate characteristic timescales, we distinguish 
between the local adjustment timescale ti and the disturbance propagation timescale t p . From Eqs. 
(PQ) and ([2]), we define the local adjustment timescale by comparing the order of magnitude of the 
terms in the induction equation. Assuming |^| ~ B/ti and |V x B\ ~ B/l, the local adjustment 
timescale is 

n = v i\i-e-'' L )- 1 , (10) 

where we have made use of Eq. In the above equation, B is the local field strength and I is 
the spatial scale over which B varies. 

The disturbance propagation timescale t p measures the time taken for a change in the photo- 
sphere (i.e. z=0) to reach a height z = h. We define it as 

where we have again assumed that |V x B\ ~ B/l. Estimates of both timescales depend on the 
choice of I, which is typically smaller at lower heights where there is a greater degree of mixed 
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polarity field. In potential field models (e.g. lGarvlll989l ). the amplitude of the Fourier coefficient at 
spatial wavenumber k scales as ~ e~ kz . So at a height z in a potential field, one may expect to 
find field gradients over spatial scales of I ~ z. If one adopts I ~ z, the local adjustment times at 
z = 1 Mm and z = 50 Mm (typical range of loop heights in an AR) would be r\ ~ 2 min and 77 ~ 7 
hr, respectively. With the same assumption / ~ z, the disturbance propagation times to heights of 
z = 1 Mm and 2 = 50 Mm are r p ~ 2 min and r p ~ 4 hr, respectively. However, these numbers for 
both the local adjustment time and disturbance propagation times should be considered as upper 
limit estimates (for the corresponding height) since we assumed that at any height, the typical 
spatial scale of the magnetic field gradient is I ~ z. While this is valid for a potential field, MF 
evolution can create sharp structures with I -C z. It is worth noting that relaxation timescales 
depend linearly on vq. For smaller (larger) values of vq, the model coronal field evolves over shorter 
(longer) timescales in response to photospheric changes and vice versa. Simulations performed with 
the present choice of vq allows the model coronal field to accumulate and store free magnetic energy 
over the course of days while still allowing for flux ropes to be created and ejected over shorter 
timescales (see section [3]) . 

The magnetic diffusivity r/ in Eq. @ has contributions from three different terms. The 
first contribution is a constant, spatially uniform resistivity of 770 = 200 km 2 s _1 . The second 
contribution has a functional form that is designed to facilitate diffusion in regions of high current 
density. It is given by 

771 " % l + exp{-(C-0.1)/0.01}' (12) 
where £ = j 2 B~ 2 A 2 and A = min{Ax, Ay, Az}. The t hird contribution is a h yperdiffusivity-like 



scheme. It is similar to the hyperdiffusivities used by ICaunt Korpil (|200ll ) in the sense that 
the diffusivity is proportional to the ratio of the third and first derivatives of the quantity being 
diffused. This type of diffusivity is efficient at suppressing oscillations at the grid scale. The form 
of the hyper diffusivity for A x is 



% y lf = ^\Al/Al\, (13) 

where A\ = j x +3 ^ 2 ^' k — 2j x + ^'^ k + j x ^^ ,k ^ A\ = j x + ^'^ k and v = m.ax.{\v\ l, ^ k , The same 

scheme is used for hyperdiffusivities in the remaining Cartesian directions. 

A second-order midpoint method is used to explicitly evolve A forward in time, with each 
update consisting of a half-step and a full-step. The adaptive time step is determined by a CFL- 
like condition given by 

At = -min{A/t; max , A 2 /??max}, (14) 

where u max and r? max are the maximum values of the magnetofrictional velocity and total magnetic 
diffusivity over the computational domain. 
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2.2. Initial condition 

For each simulation, the initial condition for A corresponds to a potential field configuration 
(i.e. current-free). This is generated by means of a potential field extrapolation using the vertical 
component of the first magnetogram in the sequence. The potential field extrapolation is performed 
with a Fourier method and assumes periodic boundary conditions in the horizontal directions. Since 
the surface field in the first frame of B z consists of mostly pre-existing weak field before an active 
region emerges, the periodicity assumption should not play a big role in determining the dynamical 
evolution of the AR as it emerges. 



2.3. Boundary conditions 

2.3.1. Top and side boundaries 

At the top and side boundaries of the Cartesian domain, the magnetofrictional velocity v values 
at the boundary cell vertices are chosen to be the same values as those defined on cell vertices one 
layer deep in the computational domain. The magnetic field B is normal at the boundaries and 
the transverse component is set to vanish. The normal component of j is symmetric across the 
boundaries and the transverse components are antisymmetric across the boundaries. 



2.3.2. Bottom boundary 

A time-dependent boundary condition based on temporal sequences of magnetograms is im- 
posed at the bottom of the computational domain (i.e. photosphere) to drive the evolution of the 
magnetic field in the model corona (z > 0). Mathematically, the choice of a boundary condition 
in this time-dependent problem involves a choice of the three components of the electric field (i.e. 
E = —dtA) at z = 0. From Eq. @, one sees that this corresponds to a choice of the photospheric 
(i.e. at z = 0) distributions of v, B, j. A consistent boundary condition, however, does not allow 
one to impose arbitrarily all three components of all three vectors. Consider, for example, the 
horizontal components of the electric field. To specify the boundary condition for these two com- 
ponents, one must specify all three components of B and v as well as the horizontal components of 
j. This choice automatically constrains the remaining component of the electric field (i.e. — dtA z ) 
since j z (z = 0) = [d x B y — d y B x ] z=0 . This examples demonstrates that an appropriate boundary 
condition for the problem permits a choice of only two out of three components of the boundary 
electric field. This property is reflected in the choice of the staggered grid. At each boundary on 
this grid, only the transverse components of A (as well as their time-derivatives) are defined. 

Having established that only E x and E y need to be specified, one still has to tackle the problem 
of how these two quantities can be retrieved from observed temporal sequences of magnetograms. 



This problem was recently analyzed in detail by iFisher et al.l 1)20101 ) . An important conclusion 
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from their work is that retrieval of the photospheric electric field from temporal sequences of vector 
magnetograms taken at one height is an ill-posed problem for which the solution is not unique. 
Where available, magnetograms at various heights in the atmosphere may help resolve the problem. 
At present however, there is no instrument capable of providing magnetograms at different heights 
at the same level of spatial coverage, resolution and time cadence that SDO/HMI can offer. Since 
SDO/HMI only provides magnetograms at one layer in the photosphere, one has to find ways to 
cope with the missing information. 

The simulations presented here are driven by electric fields derived from time sequences of 
SDO/HMI longitudinal (i.e.line-of-sight) magnetograms under a number of assumptions. Since the 
purpose of this paper is to introduce the data-driving method, the simulations presented here do not 
attempt to be fully realistic. In order to drive a simulation in a realistic fashion, photospheric electric 
fields faithful to the actual electric fields operating on the Sun must be used. The retrieval of such 
electric fields constraine d by temporal sequences of vector magnetograms (e.g. from SDO/HMI) is 



still a research problem (|Fisher et al.ll2010l . l2012bl ) and the evaluation of the quality of this inversion 



process is outside the scope of this paper. The application of electric fields derived from HMI vector 
magnetograms, and a critical assessment of the realism of the model coronal field is the topic of 
the second article in this series. 

For each active region simulation, we prepare a time sequence of input magnetograms B z , 
extracted from the series of full-disk SDO/HMI magnetograms. This quantity is used to specify 
the distribution of the vertical magnetic flux density at the bottom boundary of the Cartesian 
simulation domain throughout the simulated time period. The central position of the z = plane 
in the domain is chosen to co-rotate with the Sun at a rate (for fixed lati tude) given by the empirica l 



differential rotation profile for surface magnetic features obtained by ISnodgrass Sz Ulrichl ([1990). 
Each magnetogram in the input sequence B z spans 30.7° in heliographic latitude and longitude. The 
models here neglect the effect of curvature over the 30.7° x 30.7° patch. A plate carree projection 
is used for remapping the magnetograms. With 512 x 512 pixels in each direction, the effective 
grid spacing at the bottom boundary of the domain is Ax = Ay = 728 km, corresponding to one- 
arcsecond grid spacing in the plane of the sky (at disk center). One further assumption made is that 
the surface field on the Sun is purely radial. This assumption is used to convert the longitudinal 
magnetogram values into values for B z . To conserve flux balance in the magnetogram remapping 
process, the effects of foreshortening of surface area elements away from disk center are taken into 
account. Consecutive frames of the remapped magnetograms B z are generated at a regular cadence 
of 4.5 minutes from SDO/HMI longitudinal magnetograms. This is a lower cadence than the 90 
second cadence offered by SDO/HMI but sufficiently high to have a non- negligible impact on the 
memory and I/O requirements of the MF code. 

Given the sequence of input vertical magnetograms B z , its relation to Eh = (E x ,E y ) is given 
by the vertical component of the induction equation 
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%^ = -HVxi). (15) 

In order to solve for Eh, another relation must be specified. For example, as is done here, specifying 
the horizontal divergence of the electric field suffices: 

D(x,y) = V h -E h . (16) 

Given a chosen form of D(x,y), Fourier transforms are used to solve Eqs. (I15D and (I16p for Eh- 

Two simulations with identical initial conditions, driven by the same temporal sequence of pho- 
tospheric magnetograms B z will nevertheless yield different coronal field configurations if different 
distributions of D(x,y) are used. Since it is yet unclear how one can unambiguously and robustly 
retrieve the function D{x,y) (or Eh) from observations, the simulations in this study should be 
considered as an exploratory study of coronal field evolution under different assumptions of D{x, y). 



2.4. Method for synthetic coronal images based on a proxy emissivity 

The MF model does not treat the thermodynamic evolution of plasma in the corona. This 
means that the MF method does not provide spatial distributions of thermodynamic quantities 
such as pressure, mass density and temperature throughout the computational domain. Lacking 



age, 


Dere et al. 1997) to the synthesis of coronal images at EUV and X-ray wavelengths 


Peter et al. 


2004 




Aiouaz et al. 2005: Chen et al 


12005: Gudiksen & Nord 


und 2005: Mok et al. 


2005; 


Warren & Winebareer 


2006 




20071: Lundauist et al. 


2008a 


b: Lionello et al. 2009: 


Hansteen et al. 201C 


: Zacharias et al. 


2011 


Martinez-Svkora et al. 


2011) without making a series of assumptions. 



Instead of attempting to synthesize realistic images of the model corona using a first principles 
approach, we developed a method for calculating proxy emissivities based on the value of the current 
density-squared (J 2 ) averaged along a magnetic field line. This method allows one to conveniently 
visualize which field lines are current-carrying. 

The method works as follows. Consider a scalar emissivity field e(x, y, z) spanning the com- 
putational domain. The 'initial state' of the emissivity field is such that e = for all points. Now 
consider a point at the photospheric boundary (z = 0, i.e. bottom boundary of the computational 
domain) at position (x,y,0). Using the three-dimensional distribution of B at some time t, we 
trace this field line into the computational domain by integrating dx/B x = dy/B y = dz/B z . Define 
for each field- line trajectory the mean-squared current density 



fds, 



(17) 



where L is the length of the field line. If a field line crosses one of the side or top boundaries 
of the computational domain, we set (j 2 ) = so that the current in that field line does not 
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contribute to the emissivity. A magnetic field line will typically traverse a number of cell elements 
in the computational domain. For each of these cell elements, we increment the local value of the 
emissivity by 

de = G(j 2 )AxAy, (18) 

where G is a coefficient which can be uniform in space or be a function of B. This procedure is 
then repeated for a collection of magnetic field lines with footpoint positions equally separated by 
Ax and Ay in the horizontal directions. 

After performing the aforementioned procedure for a large number (typically of order 10 5 or 
10 6 ) of field lines, line-of-sight integrations through e(x, y, z) can be performed through any viewing 
angle to create synthetic images of coronal loops. 



3. Application to NOAA Active Region 11158 



The set of simulations presented here follow the evolution of NOAA AR 11158 over the course of 
multiple days from UTC 2011-02-10T14:00 to 2011-02-15T06:00. AR 11158 is the source region for 
the GOES X2.2 flare that took place between 2011-02-15T01:44 and 2011-02-15T02:06. There are a 
numb er of observational and theoretical studies of the X flare and its associated coronal mass ejec- 



tion (ISchrijver et al.ll201ll : Ijiang et al.ll2012l : ISun et al.ll2012l : IWang et al.ll2012l : iGopalswamy et al 



20121 ). The evolution of the AR over the days preceding the flare, however, have not received as 
much scrutiny. 

The following simulations were performed on a Cartesian grid with Ax = Ay = Az = 728 km. 
All begin with the same potential field initial condition (at UTC 2011-02-10T14:00) and are driven 
by the same sequence of magnetograms {B z ). However they differ in the choice of D(x,y), which 
we have assumed to be of the form 



D(x,y) = ttB z (x,y) 



(19) 



where ft is a parameter which is kept constant and uniform for each individual run. This functional 
form of D(x, y) is motivated by the fact that the expression for V/i • Eh = • [— v x B]h includes 
the term uj z B z , where lo z = d x v y — d y v x is the vertical component of the fluid vorticity. Thus a 
non-zero value of f2 in Eq. (|19h corresponds to imposing spatially uniform vortical motion (i.e. 
twisting motion) at the photosphere. Table [1] shows the choices of f2 for the various runs in this 
parameter study. The reference run AR1 1158170 corresponds to the case where no twisting motion 
is applied at the photosphere. 



3.1. Active region morphology 



Figure [T] shows a time sequence of synthetic coronal images and photospheric magnetograms 
(B z ) from simulation run AR1 1158170. The former were calculated using the method described in 
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Fig. 1. — Time sequence of synthetic images of coronal loops (upper panels) and underlying pho- 
tospheric field distributions (B z scaled between ±800 Mx cm -2 , lower panels) of a data-driven 
simulation of NOAA AR 11158. In this particular simulation, the twisting parameter $7 = 0. The 
synthetic coronal images (shown in logarithmic scaling) were calculated as line-of-sight integrals of 
the proxy emissivity prescription as described in section [231 
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Fig. 2. — Same as Figure [I] but for the run with twisting parameter = 1/4 turns per day. The As 
opposed to the untwisted case, synthetic images for this case show a sheared core (bright structure 
in the center of the field of view) close the to polarity inversion line between opposite polarities. 



-14- 



Run 


n 


AR11158O0 
AR111580 + 0.125 
AR11158O + 0.25 




1/8 turn per day 
1/4 turn per day 



Table 1: Properties of the various simulation runs for models of AR 11158. 17 is a parameter that 
controls the amount of uniform twisting introduced by the time evolution of the photospheric field 
(see Eq. (fI2D for definition). 

section [2~H using 1024 x 1024 field lines (i.e. four field lines per grid cell at the bottom boundary). 
The scalar field e(x, y, z) representing a proxy emissivity was calculated assuming G = constant 
(the actual value of this constant is not important for examining field morphology as long as the 
same constant is used for different snapshots, as is done throughout this article). The synthetic 
coronal images shown in the figure represent line-of-sight integrals along vertically-directed rays. 

The sequence of images in Fig. [T] gives a sense of how the model coronal magnetic field evolves 
in response to photospheric driving. At a relatively early stage of the AR emergence (left panel), 
the magnetic field lines (identified as bright loops in the synthetic coronal images) emanating from 
the AR are relatively confined. As flux continues to emerge and the AR grows, the area of the image 
covered by magnetic loops associated with the AR increases as more fieldlines become significantly 
current-carrying. 

Figure [2] shows the sequence of synthetic coronal images and corresponding photospheric mag- 
netograms for the simulation run AR11158f2 + 0.25. In this run, a twisting rate of 0, = 1/4 turns 
per day was imposed. Due to this imposed twisting, the field above the sharp polarity inversion 
line in the central portion of the AR is significantly sheared and appears as a bright core (as a 
proxy of the current content). 

The synthetic coronal images reveal that not all of the magnetic flux within the AR is internally 
connected. In snapshots for all simulation runs, one finds field lines that emanate from polarity 
patches within the AR that have conjugate footpoints in pre-existing, diffuse flux patches external 
to the AR. In fact, as the model AR grows, one finds sets of loops which seem to crawl along the 
magnetic carpet of the ambient photosphere. 

3.2. Buildup of free magnetic energy 

In the simulation run without imposed twisting (i.e. O = 0), the absence of systematic shearing 
of the AR magnetic fields results in a relatively 'quiescent' evolution of the modeled AR. This is 
reflected in the time plot of the free magnetic energy as shown in Fig [3j The free magnetic 
energy is defined as Et = E{B) — E(P), where E(B) is the energy of the magnetic field in the 
magnetofrictional model and E(P) the energy of the corresponding potential field configuration 
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(specified by the normal field through all boundaries at the same time). For the simulation run 
with $7 = 0, the free energy of the system for the duration of the simulation run is less than 7 x 10 31 
erg ~ 0.07E(P). In the run with a constant, uniform twisting rate of = 1/4 turns per day, 
Ef reaches values of ~ 8 x 10 32 erg ~ 1.2E(P). So that latter run has a much larger reservoir 
of free magnetic energy to drive ejections. This is done by means of the creation of a sheared, 
current-carrying magnetic arcade above the sharp magnetic polarity inversion line in the core of 
the modeled AR11158. Even though a series of flux rope ejections (see next section) result from 
the shearing field in the runs with imposed twisting, the free magnetic energy continues to increase. 
This is simply a result of imposing a twisting rate that is (unrealistically) spati ally and temporall y 
uniform, irrespective of the evolutionary phase of the model AR. As reported bv lJiang et al.l (|2012l ). 
the rotation patterns of the sunspots in this AR are far from uniform. 



3.3. Recurrent flux rope ejections from sheared fields 



As shown in the left panel of Fig. [H flux ropes are formed when magnetic reconnection pinches 
the upper portion of the sheared arcade. The loss of equilibrium leads to the ejection of the flux rope. 
Although magnetofriction evolves the magnetic field in a rather simplified manner, such models are 
still able to capture a number of qualitative aspects of MHD models of coronal field evolution. For 
instance, the creation of flux ropes from a sheared arcade an d their subsequent eje ction due to loss 
of forc e equilibrium is pr esent MHD simulations with shear ([Mikic fc Linker! 1 1994! ) and converging 
flows (jAmari et al.ll2000l ) about the polarity inversion line of the model bipole. F o r the particular 
case of the X2.2 fla re and the associated eruption of AR 11158, ISchriiver et al.l (|201ll ) used the 
MHD simulations of lAulanier et all (|2010l ) as an aid to interpret the observed evolution of the AR 
during the eruption. In their MHD model, the launch of the flux rope is the result of a loss of 
equilibrium when the threshold of the torus instability has been reached (jKliem Torokll2006l ). 



In the present work, the persistent twisting of the magnetic arcade in the core of the AR 
leads to the ejection of a series of flux ro pes. Similarly recurrent plasmoid/flux rope ejections were 
modeled by iManchester IV et all (2004] ), who carried out compressible MHD simulations of the 
emergence of a twisted magnetic flux tube from an idealized convection zone into a non-magnetized 
corona. The simulations showed that shear flows at the photospheric level driven by the dynamics of 
twisted emerging flux lead to recurrent flux rope ejections. In this sense, the present MF simulations 
provide a similar qualitative result. 

Another similarity between the present model and MHD models of flux rope ejection is the 
presence of a X-point type topological feature between the ejecting flux rope and the underlying 
sheared arcade. Such a feature can be seen near the location (y, z) = (0, 20) Mm in the left panel 
of Fig. [U T he combined structure consisting of the sheared arcade, X-point and flux rope were also 
reported bv lSavcheva et al.l (|2012l ). who used both an MHD model and a MF model to study the 
pre-eruptive magnetic field configuration of an AR. In their work, the formation of such a structure 
in the MHD model was the result of imposing converging flows as well as flux cancellation at the 
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Fig. 3. — Free magnetic energy for various simulation runs. Without imposed twisting (O = 0, 
solid line), the free magnetic energy is limited to ~ 7 x 10 31 erg. With imposed twisting (CI = 1/8 
turns per day and $7 = 1/4 turns per day, dashed-dotted line and dashed line, respectively), the 
free magnetic energy is up to an order of magnitude larger. 
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Fig. 4. — The imposed twisting rate of fl = 1/4 turns per day in simulation run AR11158f2 + 
0.25 leads to a sheared, current-carrying arcade which is the source region for recurring flux rope 
ejections. Left panel: Image of vertical cross-section of the normalized current distribution (|j|/|B|) 
in a vertical plane (x = —22.6 Mm) above a central portion of the simulated AR in run AR11158f2 + 
0.25. The vector field shows the transverse component of B (normalized by field strength to 
accentuate the field orientation). Right panel: Corresponding magnetogram of the vertical field at 
z = 0. [An animated version of this figure is available.] 



-17- 



polarity inversion line. In comparison, the formation of similar structures in the corresponding 
MF model was the result of the relaxation of a current-carrying flux rope inserted into an initially 
potential field configuration. 



3.4. Increase of horizontal photospheric field near the polarity inversion line 

following flux rope ejection 

An analysis of the vector magnetograms of AR 11158 taken by SDO/HMI shows a substantial 
increase in the stre ngth of the horizon tal field near the polarity inversion line after the X2.2 flare 
on Feb 15th 20 11 (IWang et alJl2012T ). In a series of NLFFF extrapolations based on HMI vector 
magnetograms, ISun et al.1 tepid )" found a downward displacement of the current-carrying channel 



above the polarity inversion line following the flare. 

As indicated by Fig.[5l a similar qualitative behavior can be found in simulation run AR11158f2+ 
0.25 after the ejection of a flux rope. In the simulation, this increase in the horizontal field strength 
results from the relaxation of the arcade field following the magnetic reconnection event which 
produced the flux rope (see also animated version of Fig. [H available online) . 

This episode of flux rope formation and arcade relaxation in the MF model occurs almost 8 
hours before the actual X2.2 flare occurred in AR 11158 so the former should not be taken as 
a faithful representation of the observed flare and eruption. The discrepancy in timing between 
the two is not surprising given the use of only line-of-sight magnetograms and the imposition of 
an ad hoc uniform twisting rate in the model. However, there is one small benefit to not having 
used the full HMI vector magnetograms to drive the simulation. This omission rules out the 
possibility that the increase of in the model is simply a side effect of the increase of B^ in 
HMI vector magnetograms. Instead, we can unequivocally attribute the B^ enhancement in the 
simulation to the relaxation of the post-eruption arcade. The same physical mechanism, namely 
photospheric response to relaxation of a post-eruption arcade (mediated by the Lorentz force, see 



Hudson et al. 



2008 



Fisher et al.ll2012al ). has been invoked bv lWang et al.l (|2012l ) to explained the 



observed increase of B^ at the polarity inversion line in AR 11158. 



4. Discussion 

In this paper, we present a framework for performing data-driven simulations of solar AR 
formation. Under the magnetofriction assumption, fluid velocities are assumed to be proportional 
to the local Lorentz force. This has the advantage that the velocity distribution is guaranteed to 
evolve the field toward a force-free state. Since the magnetic induction equation is solved to advance 
the magnetic field in time, ideal magnetofrictional relaxation preserves the topology of the magnetic 
field while decreasing magnetic energy. In regions of high current density where resistive diffusion 
(used in this Eulerian code to facilitate magnetic reconnection) is important, the magnetic topology 
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is allowed to change via reconnection. The change in the topology may permit the magnetic 
configuration to further relax to lower energy states under quasi-ideal, magnetofrictional evolution. 

By incorporating time sequences of photospheric magnetograms as boundary data, the method 
models changes in the coronal magnetic field in response to photospheric driving (including shearing 
flows and flux emergence). As an example of the application of this method to magnetogram data 
obtained by SDO /HMI, we performed a number of simulations to model NOAA AR 11158. Since the 
accurate retrieval of photo spheric electric fiel ds from temporal sequences of vector magnetograms 



is still a research problem (jFisher et al.l l2010). we choose to perform the simulations with varying 
assumptions about the underlying E that is responsible for the magnetogram evolution. In the case 
when twisting motion is absent, flux rope ejections were produced by the model. When continuous 
twisting was imposed, the sheared field above the sharp polarity inversion line in the core of the 
AR produced a series of flux rope ejections. 

Since idealized assumptions were made about the nature of the underlying photospheric electric 
field, the simulations presented in this paper are not meant to be faithful representations of actual 
eruptions from AR 11158. Nevertheless, this work demonstrates the potential utility of such a data- 
driven approach for modeling observed ARs. In a follow up paper, we will drive the model with 
actual elec tric fields retrieved f rom HMI vector magnetograms and constrained by observed Doppler 



velocities (jFisher et al.ll2012bl ). The use of electric fields more faithful to the data will facilitate a 



side- by-side comparison of the modeled AR with EUV observations from AIA (|Lemen et al.ll2012l ). 



A method for synthesis of mock 'coronal' images base on a proxy emissivity was introduced. 
The emissivity of points along a magnetic field line is assumed to proportional to the field-line 
average of j 2 . This simple technique seems to produce with a visual texture similar to EUV images 
of coronal loops. While this technique is useful for visualizing an ensemble of coronal loops in a 
simple magnetic model, it is by no means a replacement for more sophisticated techniques that use 
the thermodynamic variables from MHD models and take into account atomic physics for EUV 
image synthesis. 
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Fig. 5. — Evolution of the horizontal component of the magnetic field in the core of the modeled 
AR 11158 (simulation run with 17 = 1/4 turns per day). The field is sampled from the midplane 
of the lowest grid layer in the simulation at z = 364 km. The ejection of a flux rope between 
2011-02-14T17:55 and 2011-02-14T18:50 results in an increase of the horizontal field strength near 
the polarity inversion line. 



